Dynamics of a Tube Conveying Fluid 
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A tube conveying a large amount of fluid with a free outlet does not sit still. We construct and an- 
alyze a nonlinear evolution equation describing such phenomena. Two types of boundary conditions 
at the inlet are considered, one for which it is clamped and one for which it is hinged. Analyzing the 
linear stability of the trivial solution, we find that with the former boundary conditions, it exhibits 
a "flutter" instability, while with the latter boundary conditions, it exhibits a "rotation" instability. 
These instabilities and the nonlinear behaviors of the system are also studied numerically. 



PACS numbers: 05.45.-a, 45.30. +s, 05.45.Pq 

Thready structures appear in a wide variety of phys- 
ical systems and exhibit many different distinctive pat- 
terns and types of motion. For example, the behavior of 
chain polymers is a commonly studied topic in physics 
and chemistry. In fluid physics, the structures and dy- 
namics of vortex filaments have been studied in detail. 
Biological systems provide us with many phenomena in- 
volving thready objects from microscopic to macroscopic 
scales, e.g., the behavior of DNA, the folding of proteins, 
the beating of flagellum, and the locomotion of snakes. 
These phenomena are observed in non-equilibrium sys- 
tems in which several factors, such as elasticity, driving 
forces and dissipation, are balanced. In this class of phe- 
nomena, the motion of a tube conveying fluid with a free 
outlet has been extensively studied, owing to its simple 
nature. Paidoussis theoretically showed that the trivial, 
straight state becomes unstable at some flow rate, and 
the tube begins to flutter as the result of a Hopf bifurca- 
tion [fij. This result has been confirmed by several the- 
oretical and experimental studies |2|-||. The theoretical 
models used in these studies are physically very realistic, 
but they are applicable only for small amplitude motion. 
Also, in most studies, only one type of boundary condi- 
tions (BC) is considered, that of the cantilevered type, in 
which the tube is clamped at the inlet and free at the out- 
let. In this Letter, we adopt a rather phenomenological 
approach to construct the minimal model for motion of 
any size amplitude under more general BC at the inlet. 
We obtain an evolution law in the form of a nonlinear 
integro-differential equation. A linear stability analysis 
of the trivial solution suggests not only the existence of a 
"flutter" instability but also the existence of a "rotation" 
instability, depending on the BC at the inlet. Numeri- 
cal simulations confirm the existence of these instabilities 
and elucidate the nonlinear behavior of the equation. 

We begin by deriving the evolution equations of the 
system. For this purpose, we first clearly define the phys- 
ical system under consideration. Consider a tube convey- 
ing fluid with a free outlet. The motion is restricted to 
x-y plane, and the gravitational force is not considered. 
Also, the length of the tube is regarded as fixed. Wc 
treat the tube as one-dimensional structure and ignore 



any dependence on its radius. An elastic force acts on 
the tube in response to its being bent from a straight 
shape. In the tube, fluid flows at a constant rate, and 
its momentum change creates a force acting transversely 
on the tube. A resistive force representing some kind of 
frictional interaction with a surrounding medium is also 
included. As typical cases, we consider two types of BC 
at the inlet. One case is that in which there is a clamp at 
the inlet, keeping both the position and direction fixed. 
The other is that in which there is a hinge at the inlet, 
keeping only the position fixed. 
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FIG. 1. A discrete model of a tube. J n and S n denote 
the nth joint and segment, respectively. The Fj™ -1 ', two 
Fe i an d Fe™ +1 \ all of the elastic forces that exert on J„, 
act perpendicularly to the segments, and tend to straighten 
the tube. F„ , the force resulting from fluid flow, divides the 
angle ZJ n -i Jn Jn+i in half, and tends to increase the bending. 

We assume that a tube with the properties described 
above can be obtained as the continuum limit of the dis- 
crete articulated model depicted in Figj^. In this model, 
there are N + 1 joints and N segments. Let (x n ,y n ) be 
the coordinates of the nth joint, and let n be the an- 
gle between the nth segment and the x axis. Each joint 
has mass m(N) = M/N, where M is the total mass of 
the tube, and the segments have no mass. The fixed 
length of each segment is l(N) = L/N, where L is the 
total length of the tube. Each joint exerts an elastic 
force on itself and its neighboring joints that tends to 
straighten the tube. The magnitude of the force exerted 
is proportional to the difference between the angles of 
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the segments on both sides: = k(N) ■ \6 n +i — 9 n \ 

for 1 < n < N — 1, where k(N) is a positive coefficient 
representing the stiffness of the tube. The momentum 
change of the fluid flow at a joint creates a transverse 
force, which acts to increase the bending. The magnitude 
of this flow force is also proportional to the difference be- 
tween the angles on both sides: — w(N) ■ |0 n +i — 6 n \ 
for 1 < n < N — 1 where w(N) is a positive coeffi- 
cient representing the flow rate. Here, we assume that 
\9 n+ i — 6 n \ <C 1. The directions and the joints on which 
the forces {F^} and {F^} act are exhibited in Figj^. 
In addition to the above described forces, resistive forces 
act on each joint in the form F^ = —c(N) ■ v n for 
1 < n < N, where v n is the velocity of nth joint and 
c(N) is a positive coefficient representing the strength 
of the resistance. Since each segment has a fixed length, 
there are TV constraints. We account for these constraints 
by using Lagrange multipliers. As mentioned above, the 
inlet is either clamped or hinged. These conditions are re- 
alized by setting F^ — k9\ if it is clamped and i 7 ^ ' = 
if hinged, where F^ is a force that acts only on the first 
joint Ji. 

Our next step is to derive evolution equations for this 
discrete model. To simplify the notation, we first restrict 
ourselves to the case of a clamped tube = kOi). It 

is easy to write down 2N Newton's equations of motion 
containing N undetermined multipliers. The positions 
{xn,Hn} can be rewritten in terms of the angles {9 n } 
using the relation {x n ,y n ) — Z^™ =1 (cos6'i,sin6'j). Elim- 
inating the N multipliers from these equations yields N 
independent evolution equations for the generalized co- 
ordinates {On}. Assuming that inertial terms, such as 
those containing i: 9f, are negligible (Physically, this is 
justified for a strongly dissipative system.), we obtain the 
rather simple equation 



0(0) = d s 9(L) = 0. 



(3) 
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A ■ QijOj — B ■ Hi 



(1) 



Here, A = k/(ml), B = w/(ml), C = c/m, 9q = 0, and 
, Qij and Tlij are the components of N x N matrices 
(N-j + 1) cos(0 4 - 9j) for j > i, 



defined as Fij = Fj 
Gi 3 = Gji = (-2 + SiPf)Sij + 5 iij+1 for j > i, Hij = 
cos {(9j + 6j-i - 26»i)/2} for j > i, H io = for j < i. 

The final step in constructing our model is to take 
the continuum limit N — * oo of ([l]). In this limit, the 
discrete index n can be regarded as a continuous variable 
s G [0, L\. Carrying out an integration by parts of the 
equation so obtained, we can write the evolution equation 
for the continuous model as 



ds' / ds" {jd t 9(s")} cos {9(s) - 9{s")} 
Jo 



Here, a = lim A r^ 00 [Z 4 A(7V)], (3 = lim^oo [l 2 B(N)] and 
7 = limjv^oo C(N). These are constants independent of 
N. The convergence of these three limits is assumed, so 
that the three terms corresponding to the elasticity, the 
driving force resulting from fluid flow, and the dissipation 
are balanced. 

In the case of a hinged tube (F^ = 0) , the evolution 
equation is derived in the very same way, although the 
result differs slightly from that of a clamped tube. For 
the discrete model, in this case the above definition of Gij 
is replaced with = G ]% = (—2 + 6 iff + <5ii)<% + 6ij+i 
for j > i. Then, for the continuous model, the two BC 
are not those given by (0), but rather, 



d s 9{0) - 3 S 9{L) = 0. 



(4) 



We see below that this slight difference has a great influ- 
ence on the behavior of the whole system. 

Note that we can always set 7/a = 1 and L = 1 in 
(||) through rescaling of t and s. This means that the 
number of essential parameters governing the behavior 
of the non-inertial tube is only one, say, j3/a = e > 0. 
Note that an increase of the flow rate corresponds to an 
increase of e. This rescaling leaves (|j) as 



ds' / ds" {d t 9{s")} cos {0(a) - 9{s")} 



d*9 + esm{9{s)-9(l)}. 



(5) 



We have thus derived the evolution equations for a 
clamped tube as (||) with (||) and for a hinged tube as (||) 
with (^). In each case we set L = 1. Both systems have 
the same, trivial solution, 9(s) = 0, which corresponds 
to a completely straight tube. We now investigate situ- 
ations in which this trivial solution becomes unstable by 
considering the linearization of (^|) about this solution. 

We first consider the case of a clamped tube. Assuming 
that for all s G [0,1], \9(s) - 0| < 1 holds, we omit all 
but linear order terms in (@) to obtain 



ds' 



ds"d t 9(s") = d 2 s 9 + e{9{s) - 9{1)}. (6) 



Unlike the nonlinear equation ([5]) , the integrand depends 
only on s". Thus, we can eliminate the integrals by differ- 
entiating (||) twice. Here, in order to retain all the infor- 
mation contained in (^) , two conditions must be added to 
the differential equation that results from this procedure: 
Representing (§) as F(s) = [(l.h.s. - r.h.s.) of (§)] = 0, 
this equation is equivalent to the differential equation 
d^F(s) = 0, together with the conditions d s F(0) = 
F(l) = 0. We thus find that M) is equivalent to 



ad 2 s 9 + p sin {9(s)-9{L)}, 



(2) 



8,9 = -dt 



ed 2 s 9, 



(7) 



with two BC: 



with the two extra BC 
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<9 S 3 (9(0) + ed s 9(0) = d 2 s 6{l) = 0. 



(8) 



Now, our objective is to examine the eigenfunctions of 
(0) with the four BC given in (|J) and (|J). Expecting 
the existence of complex eigenvalues, we regard 0(s,t) 
as a complex variable 0(s,t) S C. Here, the linearized 
equation we consider is 3(0 = — 5*0 — ed 2 Q, with the 
four BC 6(0) = <9 S 3 6(0) + e«9 s 6(0) = 0.0(1) = 2 6(1) = 
0. Assuming the form 9(s,t) = e (cr+i "^$(s), with ct, w g 
.R, $ e C, we get 

_ e $" = ( CT + j w )$ 5 (9) 
$(0) = $"'(0) + e$'(0) = clamped inlet, (10) 
$'(1) = $"(l) = free outlet. (11) 

For a given e, the set of all such solutions {(a,uj, $>(s))} 
may span the solution space of (|]) with BC given in jl(i| ) 
and (|lT|). However, all we want to know at present is 
when and how the trivial solution becomes unstable. For 
this reason, we set a = in (Bh to obtain 



e$" = iw<f>. 



(12) 



We now proceed to examine the solution set of (12), 
with the BC given in ( |Io|) and (fill). First, we assume 
that the trivial solution becomes unstable beyond some 
critical value e cr . Then, there necessarily exists at least 
one solution, (e cr , ui cr , Qcr)- Note that this (e cr ,uj cr , 
need not be the critical mode itself when the a = 
eigenspace is non-simple. In fact, such case exists for 
the hinged system. 

The general solution of ( |l2| ) is easily obtained. The 
characteristic polynomial of this equation is L(P) = P 4 + 



eP 2 



0. Let us denote the four roots of L(P) 



as ±Pi(e, w) and ±P 2 (e,cu). Here, two cases must be 
considered with regard to the multiplicity of these roots. 
One case is that in which the roots are all simple. In 
this case, the general solution of ((lj) is $(s) = C\e Pls + 
C ie - Pis + C 2 e P2S + C 2 e- p2S , where Ci, Ci, C 2 and C 2 
are complex constants. The four BC given in ( |l0| ) and 
(O) for $ yield 



= 0. 



We refer to this matrix as the BC matrix B. The neces- 
sary and sufficient condition for there to be a non-trivial 
(Ci, Ci, C 2 , C 2 ) solving this equation is deti?= 0. Be- 
cause ±Pi and ±P 2 are all complex functions of e and 
ui, detB= yields two equations for e and w. There 
exists one solution, which we found numerically to be 
(e c ,oj c ) = (37.69..., 191.25...) by searching in the region 
(e,w) € [0,100] © [0,1000]. The other case is that in 
which there is a multiple root of L(P) = 0. In this 
case, the four roots are simply given by ±P± = 0, 
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±P 2 = ±i^/e = ±Pq. The general solution of ( |l2| ) is then 
$(s) =Bi + B 2 s + B 3 e p « s + B 3 e- p " s , where B 1} B 2 , B z 
and P3 are complex constants. Again considering the 
BC matrix B, we find that in this case, det£? is always 
nonzero. Thus, in this second case there is no non-trivial 
solution that satisfies the BC in ( |l0| ) and (pd|). Thus, in 
the case of a clamped tube, we have found only one solu- 
tion in which the trivial solution might become unstable. 
This is the situation of all simple roots and occurs at 
e = e c , with a Hopf bifurcation of frequency u> c . 

Next, we consider the case of a hinged tube. We apply 
an argument very similar to that above. Now we consider 
(H) with the four BC given by (0) and 

$'(0) = $"'(0) + e$'(0) = hinged inlet. (13) 

Once again, we use the BC matrix method. In this 
case, if we assume that the roots of L(P) = are all 
simple, no solution is found numerically in the region 
(e,w) £ [0,100] © [0,1000]. In the case that there is a 
multiple root of L(P) = 0, the solution set is {(e, ui = 
0, $(s) = P ) I v e > 0, V P S C}. This is the Goldstone 
mode, which exists due to the rotational symmetry of 
the hinged tube. This mode itself does not become un- 
stable, but rather is marginally stable for all e > 0. In the 
presence of the Goldstone mode, the null eigenspace may 
have a geometric multiplicity of 2. For this reason, we 
must seek a generalized null eigenfunction that satisfies 
— — e<f>" = Bq. The general solution of this equation 
is = Pi + B 2 s - {P /(2e)}s 2 + P 3 e p « s + B 3 e- p " s , 
where Bi,B 2 ,B 3 and P3 are complex constants. The 
four BC given in ( |lT| ) and (|l^) for <!> yield 
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After some calculation, we find that the necessary and 
sufficient condition for (Pi, B 2 , P3, P3) 7^ is \fi = 
tan y/e. The smallest value of e that satisfies this 
condition is Ch = 20.16..., and the null eigenfunction 
that is independent of the Goldstone mode is $(s) = 
(B /2){s 2 /e h + (cos y/eh s)/(el, cos y/eh)}. We therefore 
conclude that in the case of a hinged tube, it is possible 
that at e = the trivial solution becomes unstable and 
that a pitchfork bifurcation involving a Goldstone mode 
occurs. 

The above analysis has identified eigenmodes of the 
linearized system that may be the first to become unsta- 
ble. The results suggest that the clamped tube under- 
goes a Hopf bifurcation and begins to flutter, while the 
hinged tube undergoes a pitchfork bifurcation involving 
a Goldstone mode and begins to rotate. We confirmed 
this suggestion through the numerical simulation of (Q). 
Here, we set A = l~ 4 a, B = l~ 2 (3, and C = 7. Figure 
^| displays some typical behavior slightly above the first 
bifurcations, which are consistent with our suggestion. 
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FIG. 2. Numerical simulation of (1). (a)The clamped tube 
undergoes a Hopf bifurcation and begins to flutter. N = 10, 
a = 0.004, (3 = 0.3, 7 = 1.0. (b) The hinged tube undergoes a 
pitchfork bifurcation involving a Goldstone mode and begins 
to rotate. N = 10, a = 0.004, /3 = 0.09, 7 = 1.0. 

Figure ^ shows that the critical values e c (N), uj c (N), 
and £h(N) found numerically converge to the theoreti- 
cal values e c , uj c , and as ./V — > 00. We calculated 
e c (iV), T C (JV) = 2it/lu c (N), and e h (iV) with 1% accu- 
racy for several values of N. The data points can be 
fit to within the uncertainty on each data point by the 
function a/(N — b) + c, where a, b, and c are the fit- 
ting parameters. From this fitting we obtain the limit- 
ing values e c {N) -> 38.0 ± 0.9, lu c (N) -> 195.2 ± 4.1, 
and eh(N) — > 20.3 ± 0.6 as iV — > 00. These are almost 
identical to the theoretical values, 37.69..., 191.25..., and 
20.16.... Thus, we can say that the conjecture based on 
the results of linear stability analysis are confirmed by 
the numerical simulation. Moreover, we confirmed nu- 
merically that both bifurcations are supercritical by ex- 
amining e dependency of the magunitude of deformation. 

Let us close this Letter with some supplementary com- 
ments. If the inertial terms should be considered, there 
are two essential parameters in the equation of motion, 
rather than just one. This fact suggests that the be- 
havior of an inertial tube may be much more complex. 
Indeed, numerical simulations suggest that, while the in- 
ertial terms do not change the qualitative properties of 
the first bifurcations, they do change the subsequent bi- 
furcations. In particular, the clamped inertial tube ex- 
hibits a period doubling route to chaos, while the hinged 
inertial tube undergoes a Hopf bifurcation as a second 



bifurcation and then exhibit a period doubling route to 
chaos. Such bifurcations are not observed in the case of 
a non-inertial tube. 

An interesting feature of our model is that is relevant to 
both aspirating and discharging tubes B . In our model, 
the flow force is due only to the change in the momentum 
of the fluid at each point of the tube. The behavior ex- 
hibited by our model thus provides a prediction for both 
a discharging tube and an aspirating tube. 

The authors are grateful to Y. Kuramoto, K. Sekimoto, 
S. Toh, G. Paquette and T. Miyoshi for informative dis- 
cussions. 
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FIG. 3. The critical values of the discrete model, e c (N), 
T C (N) = 2it/ui c (N), and e h (N), plotted as function of 
l/(N - 6), where b = -0.9012, 1.9353, and 0.1126, respec- 
tively, are obtained from the fit. The intersections of the 
fitted lines and the 1 / (N — b) — axis correspond to the 
N — » 00 limit. The uncertainty of the data is much smaller 
than the size of symbols in this plot. 
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